% In this practical we will work with a single subject's data from an
% emotional faces task (data courtesy of Susie Murphy) and perform an
% time-frequency analysis in sensor space.
% This dataset can be downloaded from:
%
% www.fmrib.ox.ac.uk/~woolrich/faces_subject1_data.tar.gz
%
% Note that this contains the spm file:
% spm8_meg1.mat
% that is an SPM MEEG object that has continuous data that has already been
% SSS Maxfiltered and downsampled to 250 Hz.
%
% and
% espm8_meg1.mat
%
% which is an SPM MEEG object that has the same data epoched into the
% different task conditions.

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
%tilde='/home/mwoolrich/Desktop';
tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.2.beta.6'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

testdir=[tilde '/homedir/matlab/osl_testdata_dir'];

datadir=[testdir '/faces_subject1_data']; % directory where the data is

workingdir=[datadir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

clear spm_files spm_files_epoched;
% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[workingdir '/spmfiles/fspm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/spmfiles/efspm8_meg1.mat'];

%%%%%%%%%%%%%%%%%%%
%% DO SENSOR SPACE MULTI-BAND TIME-FREQ ANALYSIS USING OAT
% This fits the first-level GLM to sensor space data after it has been
% subject to a time-frequency transform

oat=[];
% oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[30 60]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.method='none';
oat.source_recon.dirname=[oat.source_recon.D_epoched{1} '_multitaper'];

% The hanning taper is appropriate for lower frequencies (<35Hz).  For
% higher frequencies, the dpss method allows for better control over
% frequency smoothing.
%
% http://fieldtrip.fcdonders.nl/tutorial/timefrequencyanalysis

%taper = 'hanning';
 taper = 'dpss';
 taper = 'hilbert';

oat.first_level.tf_freq_range = [30 60];
oat.first_level.tf_num_freqs = 10;

switch taper % chooses different options depending on the taper used
    case 'hanning'
        
        oat.first_level.tf_method='multitaper'; % the multitaper option allows for either hanning or dpss; both use fieldtrip's implementations
        oat.first_level.tf_time_step = 0.05;
        
        oat.first_level.tf_multitaper_taper = 'hanning';
        oat.first_level.tf_multitaper_twin  = 0.25;
        
    case 'dpss'
        
        oat.first_level.tf_method='multitaper'; % the multitaper option allows for either hanning or dpss; both use fieldtrip's implementations
        oat.first_level.tf_time_step = 0.05;
        oat.first_level.tf_multitaper_twin  = 0.25;
        
        oat.first_level.tf_multitaper_taper = 'dpss';
        oat.first_level.tf_multitaper_freqsmooth = 0.4;

    case 'hilbert'
        oat.first_level.tf_method='hilbert';

end

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

oat = osl_check_oat(oat);

%% now run the OAT

oat.first_level.name=oat.first_level.tf_method;
oat.to_do=[1 1 0];
oat = osl_run_oat(oat);

% load GLM result
stats_hilb=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

%% visualise using Fieldtrip
% note that this produces an interactive figure, with which you can:
% - draw around a set of sensors
% - click in the drawn box to produce a plot of the time series
% - on the time series plot you can draw a time window
% - and click in the window to create a topoplot averaged over that time
% window (which is itself interactive....!)

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGMAG';
S2.first_level_contrast=3;
S2.cfg.colorbar='yes';

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotTFR(S2);

%% to do a topoplot averaged over 160 to 180 ms:
cfg.xlim        = [0.16 0.18];
cfg.ylim        = 'maxmin';
cfg.zlim        = 'maxmin';
cfg.interactive = 'no';
figure; ft_topoplotTFR(cfg,data);
